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ABSTRACT 

Perturbation Theory (PT) applied to a cosmological density field with Gaussian ini- 
tial fluctuations suggests a specific hierarchy for the correlation functions when the 
variance is small. In particular quantitative predictions have been made for the mo- 
ments and the shape of the one-point probability distribution function (PDF) of the 
top-hat smoothed density. In this paper we perform a series of systematic checks of 
these predictions against N-body computations both in 2D and 3D with a wide range 
of featureless power spectra. In agreement with previous studies, we found that the 
reconstructed PDF-s work remarkably well down to very low probabilities, even when 
the variance approaches unity. Our results for 2D reproduce the features for the 3D dy- 
namics. In particular we found that the PT predictions are more accurate for spectra 
with less power on small scales. 

The nonlinear regime has been explored with various tools, PDF-s, moments and 
Void Probability Function (VPF). These studies have been done with unprecedented 
dynamical range, especially for the 2D case, allowing in particular more robust de- 
terminations of the asymptotic behavior of the VPF. We have also introduced a new 
method to determine the moments based on the factorial moments. Results using this 
method and taking into account the finite volume effects are presented. 
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1 INTRODUCTION 



Cosmology: theory - large-scale structure of the Universe - Methods: 



fluctuations below unity. These results suggest we should 
distinguish three different regimes: 



The general frame of the current theory for the large-scale 
structure formation is generally believed to be based on the 
gravitational amplification of small density perturbations. 
The detailed study of such a mechanism started in the sev- 
enties with a complete treatment of the linear growth of the 
initial fluctuations (e.g. Peebles 1980), but progress towards 
an understanding of the nonlinear aspects of the dynamics 
has been very slow. Zel'dovich (1970) proposed an inter- 
esting approximation for a qualitative description of some 
nonlinear aspects but this approximation cannot provide ac- 
curate quantitative predictions, and is certainly not valid 
in the late stage of the nonlinear dynamics. The discovery 
of the self-similar solutions with the assumption of stable 
clustering (Davis & Peebles 1977) provided a valuable in- 
sight into the fully nonlinear regime. Subsequent numerical 
analysis, however, partially confirmed the existence of such 
a regime but only for very large mean density fluctuations, 
for rms density fluctuations above ten. More recently a lot of 
theoretical and numerical efforts have been devoted to the 
study of the quasi-linear regime corresponding to density 



(i) the linear or quasi-linear regime at large scale when 
the variance is below unity; 

(ii) the intermediate regime when the variance is between 
unity and 10; 

(iii) the nonlinear regime for the smallest scales for which 
the rms density fluctuation exceeds 10. 

The nonlinear regime is expected to be reasonably well 
described by the self-similar solutions, which is very useful 
because it provides a well-defined frame for theoretical pre- 
dictions. We recall here the general results. For an initial 
power- law spectrum, P(k) oc k n , the small scale two-point 
correlation function, £2, is expected to follow a power-law 
behavior, 



6(r) 



(l) 



with an index 7 related to n (Davis & Peebles 1977). More- 
over the higher order correlation functions are also expected 
to follow a specific behavior for their global scale depen- 
dence, 



£ P (A n,...,Xr p ) = X 



-7(p-l) 



(2) 
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From this general property it is possible to show that the 
average p-point correlation functions, 



Vp 



d 3 n . 



d r p £ p (ri, 



,r P ), 



(3) 



in a cell of volume V are related to the average two-point 
correlation function with 

t P =s P %r\ (4) 

where the S p coefficients are scale independent. Note that 
£ p is identical to the p-order cumulant of the one-point den- 
sity probability distribution function. There are however no 
definitive theories for the S p parameters, although differ- 
ent models have been proposed (Hamilton 1988, Fry 1984, 
Schaeffer 1984, Balian & Schaeffer 1988, 1989a, Bernardeau 
& Schaeffer 1992, Munshi & Padmanabhan 1996). 

The intermediate regime is certainly the most poorly 
understood from a theoretical point of view. It has been the 
subject of only few semi-analytic investigations (Jain, Mo & 
White 1995, Mo & White 1996, Padmanabhan 1996, Munshi 
& Padmanabhan 1996). 

The linear and quasi-linear regime have been investi- 
gated in details over the last few years using perturbation 
theory techniques. In particular it is possible to show that 
for Gaussian initial conditions the average p-point correla- 
tion functions follow the hierarchy (Fry 1984, Goroff et al. 
1986, Bernardeau 1992), 



PT -p-i 

"p S2 i 



(•») 



coefficients (not necessarily identical to the 



where the Sf T 

S p in eq. depend on the local shape of the power spec- 
trum (Goroff et al. 1984, Bouchet et al. 1992). A lot of ef- 
fort has been devoted recently to the analytic calculation of 
these coefficients in different cases (Juszkiewicz, Bouchet & 
Colombi 1993, Juszkiewicz et al. 1995, Bernardeau 1994a). 
Thus S3 and S4 are known analytically or semi-analytically 
for a Gaussian filter and a power-law spectrum (Lokas et 
al. 1995). More interesting is the case of a top-hat filter 
for which the whole series of the coefficients S p is known 
for any cosmological model and any power spectrum in 
3D (Bernardeau 1994b) and in 2D (Bernardeau 1995). The 
fact that the complete series is known allows us to build 
the shape of the one-point density Probability Distribution 
Function (PDF). These predictions have been checked in 
peculiar cases at various levels, for the moments (Bouchet 
et al. 1992, Lokas et al. 1995, Baugh, Gaztanaga & Efs- 
tathiou 1995, Colombi et al. 1995) or for the shape of the 
PDF (Bernardeau 1994b). All these checks have been made 
for the 3D dynamics. 

Applications to observational data have, so far, given 
contradictory results. For instance Gaztanaga & Frieman 
(1994), Gaztanaga (1995) concluded that the observed 
galaxy distribution in the APM survey reproduces the the- 
oretical predictions, but Bernardeau (1995) claimed a sig- 
nificant discrepancy. In any case the fact that quantitative 
predictions, that are direct consequences of the gravitational 
instability scenario, could be checked in observational data 
boosted the theoretical investigations in this domain. In par- 
ticular the calculation of the next-to-leading order term in 
the perturbative expansion of £ p have been contemplated by 
Scoccimarro & Frieman (1995a, b), Lokas et al. (1995). In 
this paper we propose a systematic checks of the PT results 



for power law spectra, completing previous results in 3D and 
investigating the 2D case. Our objective is to circumvent the 
validity range of all the PT results at the level of the PDF-s, 
and also for the void probability function. 

The paper is divided as follows. In section 2 we recall the 
mathematical tools that originally have been developed for 
the fully nonlinear regime by Balian & Schaeffer (1989a, b), 
and turn out to be useful for the quasi-linear regime as well. 
In the subsequent section we present the specific results ob- 
tained in the quasi-linear regime, and compare them to the 
numerical results in quasi-linear and nonlinear regime. Most 
of the mathematical details about calculating S p parameters 
using method based on factorial moments are presented in 
the appendix. 



2 SCALING AND COUNTS IN CELLS 
STATISTICS 

To explore the statistical properties of the particle distribu- 
tion in an N-body simulation we consider the Count Proba- 
bility Distribution Function (CPDF), Pi(N), the probability 
of having N particles in a spherical cell of radius I. In order 
to relate the CPDF to the p-point correlation functions we 
can consider its generating function, 



(6) 



pw = E a' 

that can be shown to be given by (White 1979, Schaeffer 
1984, Balian & Schaeffer 1989a, Szapudi & Szalay 1993), 



P(\) = exp 



E 



(nV) p (X- 



Sp 



(7) 



where n is the mean number density of particles in the con- 
sidered sample. Assuming scale invariant p-point correlation 
functions we can further write, 



#(l-A)nV£ a ) 



P(A) = exp 
where the function <f> is defined by 



E 

p=i 



.S i) 



(8) 



(9) 



(we have set Si = S2 = 1). We can notice that the Void 
Probability Function (VPF) is equal to the generating func- 
tion for A = 0. As a result the probability Pi(N) can be 
derived from the VPF through the relation, 



d N P l (N) 



d\ N 



(~n) N d N P,(0) 



dn N 



(10) 



where the partial derivatives are taken at fixed volume. From 
eq. (jl^) one can show that for a continuous field 



p(S)dS = dS 



dy 



2tt£ 



exp[-0(y)/e 2 + (1 + S) y/i 2 ]. (11) 



For small £ and very small 5, p(S) reduces to the well-known 
Gaussian form, p(<5) oc exp(— 5 2 /2£), but even for moderate 
small values of 5, strong deviations from the Gaussian be- 
havior are expected (Balian & Schaeffer 1989a, Bernardeau 
1992). 
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In the nonlinear regime, under the assumption of scale- 
invariance of the correlation functions (in the sense that the 
coefficients S p are scale-independent), P(N) and p(8) ex- 
hibit characteristic scaling laws. This implies that the VPF 
has a specific scale dependence, namely that, 



a(N c ) = -]n[Pi(0)]/(nV), 

is a function of the combination 

N c = nV£ 2 



(12) 



(13) 



only. The only reasonable physical models are those for 
which a(Nc) decreases to when JV C is large (Balian & Scha- 
effer 1989a). It is thus quite reasonable to assume that, at 
large N c , a(N c ) follows a power-law behavior, 



a(N c )~aN; 



(14) 



where w is a model dependent parameter that should be 
comprised between and 1. 

Let us describe in more details the consequences of these 
assumptions in the fully nonlinear regime as they were found 
by Balian & Schaeffer (1989a). Thus, for the shape of Pi{N) 
three domains are expected. They are delimited by N = N c 
and N = N v , with 



N v = N c fa/a) 



(15) 



This number is smaller than N c when the variance is large 
(we will see that a is of the order of 1). Then, when iV < N v , 
Pi(N) is expected to follow a specific scaling, 

= (16) 
where the function g is given by 

9(z) = -z-^ f'^exp^ 1 -^ (t-t 1 -)]. (17) 

J — ioc 

It depends on to only, and when z is small it reads, 

<?(*) 



i 

z 
exp 



27RJ 





ri - w 












z 





(18) 



When N v < N < N c , Pi{N) is expected to follow a power- 
law behavior, 

-2 



P(N) = 



n c £ 2 r(w) 



(19) 



And when N > N c another scaling behavior is expected, 

Pi(jv) =^0 (2o) 

with 



h{x) = 



T^r 4>{y) exp(xy). 
2tv i 



(21) 



Furthermore, it is quite natural to expect that <j>(y) has a 
singular behavior for a negative and small value of y, 

f(y)^^s-a s T(Lj s )(y-y s )~ UJs (22) 

This singularity induces an exponential cut-off for h(x), 



h(x) 



J exp(-\y B \x), 



(23) 



when x is large, the scaling functions a(N c ), g(z) and h(x) 
have been studied in great detail computationally including 
spurious effects that may alter their measurements in finite 
iV-body catalogues (Bouchet & Hernquist 1992, Colombi et 
al. 1992, 1994, 1995). 



3 SCALING PARAMETERS FROM THE 
QUASI-LINEAR REGIME 

In the quasi-linear regime the scaling (^) for the dominant 
part of the correlation functions is a consequence of the dy- 
namical evolution. The picture obtained in the previous sec- 
tion for the shape of Pi (N) is however no more valid when £ 2 
is small. But actually the functions h(x) and g(z) can still be 
used to describe the shape of the density PDF-s. The main 
difference with the highly non-linear regime is that the low 
and large density cut-offs merge together, thus suppressing 
the domain of the power law behavior. It is then still rel- 
evant to describe the results in terms of a, to for the low 
density domain and in terms of a 3 , y s for the large density 
tail. 

We take advantage here of the fact that, for top- 
hat filtering, the whole series of the Sp T parameters is 
known (Bernardeau 1994b) through their generating func- 
tion <^>pt(j/). Note that <j>pT(y) can be seen as a low £ 2 limit 
of a more general function 4>{y,£ 2 ). More precisely it means 
that whatever y, 



Hv^) ~* <hr(y) when £ 2 -* 0, 



(24) 



but there is no guarantee that this convergence is uni- 
form, that is that the limiting function <j>pT(y) is reached 
at the same time in terms of low values of £ 2 whatever y. 
It is however strongly suggested by the numerical results 
(Bernardeau 1994b, Colombi et al. 1996), and we will make 
this assumption in the following, thus considering the global 
properties of </>pT(y) as a valid model for the calculation of 
the scaling functions, a(N c ), h(x), g(z). 

So the function, <?W(y) = ~Yl™=i S p T (-J/) p /p ! - is 
given by the system, 



4>PT{y) 

r{y) 



y + yGpT[T{y)] + -T (y) 



-y 



dg PT [-r(y)] 
dr 



(25) 



where the function <5pt(t) can be deduced from the spheri- 
cal model dynamics, since we have 



<5pt(t ) = Gs 



<t m (Mo[1+0pt(t)]) 



ctm(Mo) 



(26) 



where Gs (t) gives the quasi-linear density contrast as a func- 
tion of the linear density contrast, r, and om is the rms den- 
sity fluctuation at a given moss scale, Mo being the mass 
scale associated with the filtering radius. These results are 
valid for the 2D and 3D dynamics. For power-law spectra, 



P(k) ~ k n , 

the generating function Gpt(t) reads, 
Gpt(t) = Gs (t [1 + GpT(r)]-^ +d)/i2d) ) 



(27) 



(28) 



To get quantitative predictions from PT one needs to know 
the expression of Gs{t) in both cases. There are no simple 
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analytical expressions for Qs(t). It is however possible to 
develop Qsir) with respect to r and one obtains, 



r r \ j- 17 2 341 3 

Qs(r) = -r+-r -— r + ... 

for the 3D case (Bernardeau 1994b) and 
29 3 



Ss(r) 



, 18 2 

-T H T 

21 



42 



(29) 



(30) 



for the 2D case (Bernardeau 1995). One can then compute 



si 



s: 



3D 



34 _ 
7 

6712 



(n + 3), 
62 



1323 3 
and for the 2D case, 



(n + 3) + J.(n + 3) 2 , 



C 2D 
^3 



36 3 , 

y-2 (n + 2) - 

2540 



49 



-33(n + 2) 



21 I 

4 V 



2) 



(31) 
(32) 

(33) 
(34) 



To have more global properties of (j>(y) it is necessary to 
know the global shape of Qs{t) and one can actually show 
that, 



fc(T)=(l-l)"" 

with 

3D, 

2 
and 

13-1 



- 1, 



v = — for 
2 



1.3 for 2D, 



(35) 



(36) 



(37) 



is actually a good approximate function. The parameters 
v are chosen to reproduce the exact asymptotic behavior of 
Gs(t) + 1 for large r. Note that the Zel'dovich approximation 
would have given v = d in (^), (Bernardeau & Kofman, 
1995). It is interesting to remark also that the asymptotic 
behavior of 4>(y) for large values of y, and thus the value 
of u, is given by the asymptotic behavior of Qs(t) + 1 for 
large r (Bernardeau & Schaeffer 1992). From (j28h it is easy 
to show that 



Gpt(t) + 1 
so that 



_-v / (l-v(n+d) / (2d)) 



(I 



d(2 + v)/u-(n + d) 

As a result for the 3D case we have, 

3D 3 

7- (n + 3)' 

and for the 2D case, 



2(3 + VT3)/(V13 - 1) - (n + 2) 



(38) 



(39) 



(40) 



(41) 



5.1 -(ra + 2)' 



Note that these results are exact, in the sense that they 
do not depend on the approximation (^B|) made for Qs{t). 
These parameters entirely determine the shape of the func- 
tion g(z). It can be however more interesting to consider a 



Table 1. Parameters of the singularity, eq. (|43|), for the 2D PT 

case 



11 


Vs 


<ps 




a 3 


-2 


-0.171979 


-0.19731 


1.60005 


-1.71747 


-1.0 


-0.211979 


-0.251915 


l.oUDUZ 


-Z.Z4olo 


-1 


-0.276897 


-0.349585 


2.22626 


-3.40954 


-0.5 


-0.402865 


-0.580901 


3.54813 


-7.73243 


Table 


2. Parameters of the sin 


gularity, 


J — i 

eq. (^3|), for 


case 










n 


Vs 


<t>s 






_3 


-0.1848 


-0.214286 


1.6565 


-1.83824 


-2.5 


-0.213447 


-0.253791 


1.80389 


-2.20845 


-2 


-0.252692 


-0.3113 


2.0344 


-2.80188 


-1.5 


-0.309854 


-0.402947 


2.44269 


-3.93026 


-1 


-0.401244 


-0.572949 


3.3437 


-6.66931 


-0.5 


-0.57352 


-1.00758 


6.63172 


-18.6473 



more general expression given by the saddle point approxi- 
mation in eq. dll 



p(S)dS 



d5 



1 - rGp T (T)/gp T (T) 
2na 2 



1/2 



exp 



2<j2 i ' with ^ PT ( r ) 



(42) 



This expression reduces to (jl^) when the functions Q PT (r) 
and Gpt(t) axe replaced by their asymptotic power-law be- 
havior. 

To obtain the position of the cut-off in h(x) one needs 
to know the dominant singular value of <^pt(j/) and the be- 
havior of 4>pt(jj) near this value. It is actually not given 
by the singularity appearing in the expression of Gs(t) but 
generically by the 2 nd equation of the system (|25|). It is quite 
easy to see that there is a value of y for which (dy/dr) = 0. 
At this point we have 

\3/2 



<Vr(y) = 4>s + r a (y - Vs) + a s (y- y s }"'" + 
As a result we get an asymptotic shape for h(x) 



h(x) 



3a s a 



(S + l-r 3 ) 



-5/2 



4v^F 

exp[-|y s |(l + «5)/a 2 + |^|/a 2 ]. 



(43) 
2;iven by , 

(44) 



The parameters <j> s , a 3 , y a and r a are given in table 1 for the 
3D case as a function of the power law index n and in table 
2 for the 2D case. 

One can see that the singularity is sharper for low values 
of n, and for the 2D case. Note that for n > there is no 
singularity anymore, and the form (^) only can be used to 
describe Pi(N). In such a case the asymptotic behavior of 
Gpt{t) is 



2d/(n+d) 



(45) 



1 + Gpt{t) w — 

so that \t c \ ~ 1.47 for the 2D case and |t c | ~ 1.69 for the 
3D case. As a result the large density tail takes the form, 

t c (n + d) 



p(<5)d<5 



2d 



(1 + 5) 



(d-n)/(2d) 



(46) 



* In Bernardeau (1992), the regular term. r s (y — y s ), was ne- 
glected so there is a slight change in eq. (Uj) 
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Xy / 2n/(n + d) exp(-^(l + <5) (n+d)/d ) dS 
which gives a sharper cut-ofF than in the expression (Q) . 



4 COUNT PROBABILITY DISTRIBUTION 

FUNCTION IN NUMERICAL SIMULATIONS 

4.1 The simulations 

The simulations used here are numerical models for the grav- 
itational dynamics of collisionless particles in an expanding 
background. We are studying evolution of initial Gaussian 
perturbations in fi = 1 universe. All the simulations are 
done with a particle-mesh (PM) code with 512 2 particles 
in 2D in an equal number of grid points and 128 3 particles 
with 128 3 grid points in 3D (Melott 1986; Melott, Weinberg 
& Gott 1988, hereafter MWG). More details about the pe- 
culiarities of the simulations used here can be seen in Melott 
& Shandarin (1993). The code has at least twice the dynam- 
ical resolution of any other PM code with which it has been 
compared. 

We use a subset of initial conditions used by Beacom et 
al. (1991) and Kofman et al. (1992) in their studies. The time 
evolution of the N-body models can be seen in the video ac- 
companying the paper of Kauffmann & Melott (1992). The 
models for which we carry out our analysis are featureless 
power-law spectra of the general form, 



P(k) oc k n for k < k c , 
= for k > k c . 



(47) 
(48) 



We have analyzed power-law models with n —2, 0, -2 in 2D 
and n =1, 0, -1, -2 in 3D with a cutoff in each case at the 
Nyquist wave number: k c — 256 kf for 2D and k c = 64 kf for 
3D where kf = 27r/Lbox is the fundamental mode associated 
with the box size. 

We choose <t(A;nl), the epoch when the scale 27r/feNL is 
going nonlinear as a measure of time. 



J^ y P{k)kdk' 

Cr(feNL) = I — r 

' J^P(k)kdk 



(49) 



The first scale to go nonlinear is the one corresponding 
to the Nyquist wave number. This happens, by definition, 
when the variance a is unity. Of course as a increases succes- 
sive larger scales enter in the nonlinear regime. The simula- 
tions were stopped at Anl = 2/g r id, 4i gr id, 8£grid, Lbox/2. 
The growth rate of various modes in linear theory were stud- 
ied in MWG for this PM code. The results given by our code 
at A = 3Z g rid are equivalent to the ones obtain by a usual PM 
code at A = 8£ gr id- This is due to the staggered mesh scheme. 
So we expect that our code performs well at the wavelength 
associated with four cells and since the collapse of 4 /grid- 
size perturbations will give rise to condensations of diameter 
2 Zgrid or less, the smallest cell size that we take into account 
should be bigger than 2 Z gr id ■ On the other hand Kauffman 
and Melott (1992) found that for voids of size greater than 
size Lbox/4 self-similarity was broken in a model equivalent 
to our index n — — 1 in 3D, see also Gramann (1992) and 
Melott and Shandarin (1993). We therefore restrict our cell 
sizes to be less than L/10. As a result our cell sizes vary 
between 2L gri d and Lbox/10. We also do not use cells with 



a < 0.1, since in this case shot noise becomes comparable 
to the fluctuation power impressed on the simulation. 

Our simulations were started by using the Zel'dovich 
approximation (Klypin & Shandarin 1983) but we wait 
long enough before doing a comparison with the PT re- 
sults so that the effects of the Zel'dovich approximation have 
died away (see for instance Baugh Gaztanaga & Efstathiou, 
1995). 



4.2 The Count Probability Distribution Functions 

To evaluate the CPDF we calculate the occupancy of spher- 
ical cells of size I disposed on a regular mesh. The number 
of cells to be used is to be as large as possible so that all 
structures are fully taken into account. This can be achieved 
by considering cells that are about l c apart where l c is the 
typical separation between particles in clusters (N c (l c ) — 1). 
Actually the major constraint comes from the resolution of 
the iV-body code. If l les is the resolution scale of the N-body, 
it is clear that probabilities below (l ICB /L) d are meaningless. 
As a result we probe probabilities as low as 10 -6 both in 2D 
and 3D. 

Finite size of the sample affects mainly the large N tail 
of the CPDF, and this effect is all the more important that 
the cell size I is large. The main effect is that the large 
density tail of the PDF is dominated by a single cluster. It 
creates a bump in the PDF which is followed by an abrupt 
cut-off of the distribution. These features have been rec- 
ognized as faked by Colombi et al (1995) and methods to 
correct for it have been proposed. This is important in par- 
ticular for the derivation of the high order moments that are 
extremely sensitive to the defects in the large density tails. 
In the method proposed by Colombi et al., this effect is cor- 
rected using a theoretical prejudice, that is that the density 
PDF is assumed to follow an exponential cut-off. This hy- 
pothesis is supported by the PT theoretical results and by 
our understanding of the nonlinear regime. We will use this 
method as well to compute the high order moments. 

We have followed the evolution of CPDF for all the 
spectra for different scales with time. The different scales 
that we have studied are in the range —2.2 < log(//Lb ox ) < 
— 1.2 and they are separated by equal logarithmic intervals 
Alog(Z/L b ox) =0.2. 



4-2.1 Results in the Quasi-linear Regime 

In figs. 1 and 2 we present the results of the PDF-s for the 
2D and 3D dynamics for different values of the variance. 
The latter is at most of the order of 1.5 (see table) and the 
agreement is found to be extremely good for the smaller 
variances. For variances approaching unity, the departure 
from the PT results depends on the value of the initial index. 

When there is a lot of power at large scale (n small), 
the agreement is better then for the other case. This is not 
too surprising since when n is large there is a lot of power 
at small scales that can affect the behavior of the largest 
scales. 
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Table 3. Parameter values for the P(N) in 2D dynamics 



n = 2 


n = 


n = -2 




f 0.61 


0.47 


0.88 




£ 1.04 


1.00 


1.32 




£ 1.58 


1.39 


1.82 




nv 4.11 


10.33 


69.55 




nv 1.63 


1.63 


4.11 




nv 0.65 


0.65 


0.65 




Table 4. Parameter values for the P(N) in 


n = 1 


n = 


n = -1 


n = -2 


£ 0.40 


0.48 


0.57 


0.71 


£ 0.79 


0.86 


0.97 


1.09 


£ 1.44 


1.48 


1.5 


1.64 


nv 33.23 


33.23 


33.23 


33.23 


nv 8.33 


8.33 


8.33 


8.33 


n?) 2.09 


2.09 


2.09 


2.09 



4-2.2 Results in the nonlinear Regime 

The nonlinear regime has been explore mainly in terms of 
the function h(x) (see eq. Ji^]). Indeed when the variance is 



large it is natural to expect that N c £ 2 Pi (N) is a function 
of N/N c only (for a given power law index). As a result one 
expects that, when plotted with the appropriate variables, 
the PDF measured at for different smoothing scales coincide. 
This test is presented on figs. 3 and 4. Here we see that 
the locations of the PDF are indeed the same when the 
adequate change of variable is made. It should be noted 
however, that not all the PDF have been plotted. As the 
function h is pertinent is the large density tail only, the 
PDF-s have been truncated in the law density domain. We 
have removed part of P(N) where it is dominated by shot 
noise. Typically N < 10 while increasing with large scale 
power. Very large TV part of P(N) is dominated by large 
statistical fluctuations, we applied smoothing to reduce such 
fluctuations. 

In order to have quantitative results we used the param- 
eterized fit proposed by Bouchet, Schaeffer & Davis (1991) 
to describe the function h(x), 



h(x) 



a(l 



! exp(-x\y s \) 



r( w ) 



(1 + bx) c 



(50) 



The values of u) and a are estimated from CPDF (see in 
the next subsection), \y 3 \ is found from fitting the large iV 
tail of CPDF which shows an exponential cut-off. The other 

© 1996 RAS, MNRAS 000, 000,000 



bcaling in Gravitational Clustering 7 





1 
i 











r- 


C\J 




X 


-1 


o 




m 


-2 








-3 




n = 2 i 



-1 1 

!gio x 



x 



1 



-1 



00 -2 
-3 



n = \ 



-1 o 



■10 



1 

x 



X 



1 



-1 



tail -2 
-3 




n = — 2 



-1 o 1 

Igio x 



Figure 3. The measured h(x) in 2D N-body simulations for n = 2, n = and n 
correspond to fitting function of the form (50) with the parameters of table 5. 
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Figure 4. The measured in 3D N-body simulations for n = 1, n = ,n = —1 and n = —2 in highly nonlinear regime. Long dashed 
lines correspond to fitting function of the form (50) with the parameters of table 6. 



Table 5. Parameters of the fitting function h(x) (eqj50|) for the 

2D case 



n 


U! 


a 




b 


c 


2 


.85 


3.35 


1.31 


3.2 


-.83 





.72 


2.51 


0.58 


.06 


1.2 


-2 


.30 


1.2 


0.09 


1.1 


.80 



Table 6. Parameters of the fitting function h(x) (eqj50[) for the 

3D case 



n 




a 


\Vs\ 


6 


c 


1 


.70 


2.38 


.64 


.00 


.00 





.55 


1.59 


.11 


.38 


.60 


-1 


.40 


1.31 


.23 


.65 


.70 


-2 


.33 


1.25 


.11 


.80 


.95 



two parameters are adjusted to reproduce the constraints 
Si = 1 and & = 1 to at least 5 percent accuracy. 

We give the parameter values that have been used in 
tables (g |). 

The variation with the power law index is extremely 
large, in particular for the value of uo. It is interesting to 
note that the results are better when n is large, for which 
there is a convincing overlapping of curves. The situation 



is more questioning when n is small, but it is still not clear 
whether it is due to some numerical difficulties (finite volume 
effects are large when the power law index is small) or to a 
genuine physical effect. In particular it has been found that 
the S p parameters reach their asymptotic value for quite 
large values of the variance. The curve presented here may 
therefore still be in the intermediate regime for which the 
function h(x) is rapidly changing. 

We can also note that the CDM case as analyzed by 
Bouchet, Schaeffer & Davis (1991) matches roughly with 
our calculations ofn = — 1 or n = - 2 for 3D dynamics with 
their computed value of u — .4 ± .05 and \y s \ = 0.08. 



4.3 The Moments 

Once CPDF has been calculated one can use this informa- 
tion to compute moments of this distribution and hence S p 
parameters (eq.[4]). As recalled in the introduction these 
parameters have been studied extensively in recent past, 
both analytically and computationally. In the limit a 2 — *• 
and for Gaussian conditions they are constant (Bernardeau 
1992) and can be calculated using perturbation theory. The 
resulting values of 5*3 and 54 are given in the first section. 
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Figure 5. The measured S p parameters in 3D N-body simulations for n = 1, n = 0, n = —1 and n = —2 are plotted against £2- The 
filled dots our measurements of S p parameters for p = 3 to p = 5 after taking all the corrections. Crosses represent results from Lucchin 
et al. (1994). Open circles correspond to measurements of S p parameters by Colombi et al. (1996). Solid lines represent predictions from 
perturbation theory. 



Figure 6. The measured 56 parameter after (filled circles) and before (open circles) finite volume corrections for different initial power 
spectra n. 
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Figure 7. The computed S p parameters are plotted against £2 for n 
predictions from perturbation theory. 
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Numerical results show good agreement with theoretical pre- 
dictions. 

Of course these formula are not expected to be valid in 
the intermediate and in the nonlinear regime. Note however, 
that the spherical model might be of help to get insights into 
those two regimes as pointed out by Mo & White (1996), 
Munshi & Padmanabhan (1996). 

In this paper we have developed a new method based 
on factorial moments (see Appendix) to compute S p pa- 
rameters, with corrections due to finite volume effects taken 
into account. Using this method we have computed S p pa- 
rameters up to p — 6 for both 2D and 3D and compared 
them with earlier results for power law spectra. We found 
that finite volume corrections are extremely important as 
scale of nonlinearity increases. They are most important for 
spectra with lot of large scale power e.g. n — —2. Com- 
parison of our results with earlier results of Colombi et al. 
(1996) where different method of volume corrections were 
applied shows reasonable agreement within the dynamical 
range studied by us. We however find that the S p coefficients 
in the non- linear regime are reasonably constant: the drift 
S p oc £ 0,045 ' p 2 seen by Colombi et al. (1996) is no longer 
present. The difference in the finite volume corrections ex- 
plain why the agreement is better for spectra with less large 
scale power. To do these corrections we are a using a sys- 
tematic method - see Appendix - which makes sure that all 
normalization constraints are well respected. Note however 
that the corrections depend mainly on proper parameteri- 
zation of h(x). With increasing large scale power accurate 
determination of h{x) and a(N c ) becomes difficult. This may 
also partially explain the lack of agreement between our re- 
sult and Colombi et al. (1996) for n = — 2 spectra. 

The S p parameters were also studied by Lucchin et al. 
(1992) using central moments. Our results shows good agree- 
ment for all spectra with their results. Although finite vol- 
ume corrections were not taken into account they used sev- 
eral realizations of the same spectra which made their results 
agree with ours where volume corrections were taken care of 
and only one realization has been used. Since variation of S p 
was not studied for large values of £2 in Lucchin et al. (1994) 
it is difficult to compare their results with our results. 



4.4 The Void Probability Function 

In the previous section we have investigated mainly the 
large density tails. As mentioned earlier the void probability 
function is directly related to the clustering properties but 
contains complementary information in numerical measure- 
ments. The scaling in the CPDF statistics is related to the 
scaling in a = — hi[P(0)]/(nV). Note that in the absence 
of clustering, for a Poisson distribution, a is exactly unity 
and when some clustering is present the VPF is expected to 
grow (it is more probable to find a void) and hence a de- 
creases. Moreover from the previous sections we know that 
the function a when expressed as a function of N c is ex- 
pected to have a power-law behavior for large iV c . For small 
N c it tends to the Poisson limit, a(N c ) — 1. 

We have considered seven different epochs starting from 
the epoch when the grid scale goes nonlinear to the epoch 
when the box scale goes nonlinear. We have also considered 
different level of dilution to the full set of 512 2 and 128 3 
data. 

Different spurious effects affect the VPF differently, un- 
like CPDF it is less affected by finite volume of the sam- 
ple but it is very much affected by grid effects. Over-dense 
regions lose memory of the initial grid just after the first 
shell crossing but in under-dense regions the grid structure 
is present till very late stages and these under-dense regions 
contribute more to VPF. It was shown that this effect is 
negligible if we restrict ourselves to those cells where the 
conditions P(0) > 1/e is satisfied (Colombi et al. 1995). 
Grid effects are more significant for spectra with more power 
on larger scales where as for spectra with lots of power in 
smaller scales collapse of smaller objects erases the memory 
of initial grid very fast. 

Scatter in plots increases from n = 2 to n — and for 
n — —2 the a(N c ) diagrams reveal an significant evolution 
with time. For n = — 2 spectra larger and larger modes have 
more power and results are affected badly by finite volume 
corrections. Since we superpose all scales from quasi-linear 
to highly nonlinear regime in the same plots, any scatter 
gives a measure of variation of S p parameters in these two 
regimes. In earlier studies, it was noticed that variation of 
S p parameter from quasi-linear to highly nonlinear regime 
increased with presence of large scale power. This may ex- 
plain the small scatter in our plot for n— 2 spectra. We do 
not see such a trend in our 3D analysis which may be due to 
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Figure 8. The measured <j(N c ) in 2D N-body simulations for n = 2, n = and n = —2. The open symbols correspond to cases where 
the variance is below unity, filled symbols for a variance above unity. The slope of the solid lines are given in table 7. 



Figure 9. The measured cr(N c ) in 3D N-body simulations for n = 1, n = 0, n = —1 and n = —2. The open symbols correspond to cases 
where the variance is below unity, filled symbols for a variance above unity. The slope of the solid lines are given in table 8. 



Table 7. Measured values of u> for the 2D dynamics, Fig. 6, 
compared to the PT prediction, eq. (39). 

n = 2 n = n = —2 
^moas. 0g5 0.72 0.3 

w PT - 0.65 0.39 

Table 8. Measured values of u for the 3D dynamics, Fig. 7, 
compared to the PT prediction, eq. (38). 

n = 1 n = —1 n = n = 1 
^moas. Q 7 55 04 033 

w PT 1 0.75 0.6 0.5 

the smaller available dynamic range. Computation of a(N c ) 
for very small values of £2 is constrained by the restriction 
P(0) > 1/e which is to be satisfied for avoiding grid effect 
and discreteness effect. 

The values of 10 we get are given in table (Q) for the 2D 
case and in table Q for the 3D case where the errors on the 
measured values are about 0.05. It is interesting to see that 
PT predicts the right trend for the n dependence, although 
there are significant discrepancies between the measured u 
and the ones predicted by PT. 



5 CONCLUSION 

We have shown that PDF constructed from PT formalism 
works extremely well for a 2 < 1. Since all the loop correc- 
tions to S p parameters are neglected in this kind of approach 
it seems that loop corrections may not be very important in 
the perturbative regime for construction of PDF. 

The constructed PDF is more accurate for spectra with 
less power on smaller scales. The existence of lots of small 
scale power produces a flow of power from smaller scale to 
larger scales which contradicts the basic assumptions of per- 
turbation theory where density evolution at sufficiently large 
scales can always be described by linear theory. 

For spectra with more power in larger scales the evolu- 
tion is quite rapid and Pjv take the characteristic nonlinear 
power- law form quite early even when a < 1. 

We developed a method based on factorial moments to 
calculate higher order correlation functions and used them to 
study evolution of S p parameters for power law spectra in 2D 
and 3D. Comparison with earlier studies shows reasonable 
agreement. 
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APPENDIX 

5.1 P(N) and its continuous analogue H(v) 

Using the relation between count in cell and void probability 
function one can write, 



P(N) 



d 



■ exp 



N\ d/j, N ~" c \ & 
which can also be written in the following form 

0((1-A)7V C )' 



P(N) 



1 



~\n+T exp 



6 



(51) 



(52) 



The above integral has to be evaluated along a contour 
around A = 0. 

We then define the function ^(t) — —cj>(—t) with which 
the function is defined, 



exp 



1st 

^exp(— )IL{v). 



(53) 



' £2 Jo 

Using the definition of P(N) now one can easily show that 



P(N) 



(54) 



Therefore can be viewed as the continuous limit of 
P(N) in the limit of large number densities. This can be seen 
by change of variable A = 1 + t/N c in equation (2) and then 
taking the limit N c — > 00 and N — > 00 with the ratio N/N c 
remaining finite which gives P(N) = n(A^). More precisely 
P(N) is the convolution of the function with a Poisson 
distribution describing the shot noise effects. 



5.2 Factorial moments and S p parameters 

The factorial moments of P(N) can be related with moments 
of by the following expression. 



^2 N ( N - -p + l)P(N) 



v v Yl(v)du. (55) 



Now let us expand 9(t) in a power series of t, ^(t) = 
y"]^ \& p t p . One can convince himself that the coefficients ^ p 
have the following relation with the S p parameters; S p = 
p\^ P - Since in realistic scenarios, S p behaves as p\ the 
coefficient are expected to be of order unity. 

Similarly one can define the normalized factorial mo- 
ments of P{N) by the following expression 



6 



' N(N - 1)...(JV -p+ l)P(N). 



(56) 



These numbers are also of order unity. One can define a gen- 
crating function E(t) for E p parameters E(t) = ~E p t p . 
We have the following expression connecting these two gen- 
erating functions ^(t) and E(f), 
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! + ^ =exp(-^), 
which can also be written as, 
*(t) = f 2 ln(l + H^). 

?2 



(57) 



(58) 



Expanding the above relation in powers of l/£ 2 we can write, 

(59) 



T . . ^, . lE 2 (i) lS 3 (i) 



2 6 3 ei •• 

One can then write down the expressions of \]/ p 's as a 

function of E p 's. We present some lower order relations here, 

*i = Ei (60) 

*2 = E 2 -iSi (61) 

£ ?2 

»^E 4 -(E,E, + l S ^ + 5|2-i| ,63, 



* 5 = E 5 - (E1E4 + E 2 E 3 )-^ 

+ (E?E 3 + E 1 E 2 )|-^ + l| (64) 



* 6 = E 6 -(EiE 5 + E 2 E 4 iE^)i- 

2 42 



+ (E?E4 + 2E 1 E 2 E3 + is^)^ 

3 z\ 



These formulae are valid for arbitrary values of £ 2 . Note that 
of course in the limit £ — > 00 we have E p = S p . 



there are (l/l c ) 3 grids and (L/l c ) 3 cells, but only (L/l) 3 cells 
are completely independent. 

Events corresponding to P(N)'s smaller than (l/L) 3 
will be either over represented in case there is one such 
event in the sample, or under represented in case there is 
none in the sample. So it is clear that one has limited access 
to P(N)'s smaller than the above value. From the exponen- 
tial decay of TL(v) oc exp(— v/v*) with z/* — x,N c one gets a 
limit for TV above which the information about P(N) is not 
contained in the sample. 



N max w x,N c \n(L/l c f 



(66) 



A more realistic parameterization of II(i/) for v <C TV C would 
be 



litis) oc —Lfo(Jl) 
which leads to the formula 



TV„ 



x*N c In 



L\ yfrx* 5/2 / L 



(67) 



(68) 



For JV > N m ax, P(N) abruptly drops to zero due to the 
finiteness of the sample. The moments E p are dominated by 
values of v that can be inferred to be 



^ ~ (p- 5/2)i/,; 



p< 5/2 



(69) 



These coefficients are not known for z/e > N max that is 
for the simple for (19) p must satisfy 

p < 5/2 + 3 ln(L/Z c ) (70) 

P(N) whose value is between (l/L) 3 and (l c /L) 3 con- 
tains systematic wiggles due to the fact that only (L/l) 3 
cells out of (L/L c ) 3 cells are fully independent but averages 
of P(N)'s such as the moments (5) are less sensitive to this. 
It is nevertheless better to use this information than to drop 
it. With (L/l) 3 trials P(N) will become inaccurate much ear- 
lier. E p are systematically under estimated for large p due 
to abrupt drop in P(N) for JV > N max . 



5.3 Minimum and maximum length scales to be 
probed 

The particle positions in simulations sample an underlying 
continuous field. To extract this field from actual output, 
one has to work with sufficiently large cells, so that the field 
points appear continuous. The condition is given by N c ^> 1. 
In practice if l c is the scale where iV c = 1, scales of few times 
l c start to be usable and obviously the code resolution pro- 
vides the minimum length scale that can be probed. Choice 
of maximum length scale is slightly arbitrary as larger and 
larger cells starts finite volume effects. These effects can be 
visualize however with irregularities in the shape of the mea- 
sured PDF-s. 

5.4 Number of trials 

The typical distances between particles in a cluster is l c . If 
the sample is divided in cells of size I some information at 
scales smaller than I is erased. To recover all informations 
available in the sample one has to use many grids displaced 
by a distance l c from each other. As mentioned earlier in 3D 



5.5 Correction for finite volume effect 

A simple way to correct E p is to use the form (17) to supple- 
ment the missing information at large JV, and to use P(N) 
for TV < N max . 

P corr (N) = (l-^-b {N ~ Nc) )P com (N); N<N max 

?2 "c 

P corr (N) « n(z/); N>N max (71) 

It is clear that the corrected P(N) so constructed have 
to be normalized properly before using it to calculate S p 
parameters. 

Using the constrain ^ P(N) = 1 we get a = Ho and 
~^2NP(N) = 1 gives us b = Hi. Where we have used the 
following notation 

/•oo 

H p = / x p h(x)dx (72) 

J N maal /N c 

Where we have neglected the corrections of order 1/^2- 
The corrected 7V C now can be written as 

N corr = { i_ QHl Sg + H 2 )K° mp (73) 
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and finally the corrected E p 's are of the form 
E™ p - (p + 1)H 1 Y; C °™ P + J 1 H„ 

ycorr _P v+ P' (7d\ 

p ~ (l-6ffiE 3 +ff 2 )i>-i ( ' 

Corrected S p parameters can now be recovered by us- 
ing relation between ty p and E p as described already. This 
way to correct has the advantage of being rather simple and 
preserve all normalizing conditions. 

This method of calculating S p parameter is an alterna- 
tive to the method generally used based on central moments. 
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